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This paper discusses the application of the constraint force equation methodology and its implementation for 
multibody separation problems using three specially designed test cases. The first test case involves two rigid bodies 
connected by a fixed joint, the second case involves two rigid bodies connected with a universal joint, and the third test 
case is that of Mach 7 separation of the X-43A vehicle. For the first two cases, the solutions obtained using the 
constraint force equation method compare well with those obtained using industry- standard benchmark codes. For 
the X-43A case, the constraint force equation solutions show reasonable agreement with the flight-test data. Use of the 
constraint force equation method facilitates the analysis of stage separation in end-to-end simulations of launch 
vehicle trajectories. 
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Nomenclature 

rigid body A and rigid body B 
joint location in body A and body B 
components of sensed mass center 
acceleration in local body frame, ft/s 2 
joint constraint force vector for body A and 
body B, lby 

external force vector acting on body A and 
body B, Vctf 

constraint force components in local body 
frame, Vo f 

external force components in local body 
frame, Vo f 

inertia tensor about mass center for body A 
and body B, slug • ft 2 
mass of body A and body B, slug 
number of equations and unknowns required 
to solve constraint force equation problem 
inertial position vector of A and B, ft 
joint constraint torque vector for body A and 
body B, ft - lby 

external torque vector acting on body A and 
body B, ft • lb f 

constraint torque components in local body 
frame, ft - lby 

inertial components of mass center velocity 
relative to inertial reference frame, ft/s 
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*A> x fl 

= inertial position vector to mass center of body 
A and body B, ft 

*A. x fl 

= linear acceleration vector of body A and body 
B mass center, ft/ s 2 

a 

= angle of attack, deg 

P 

= angle of sideslip, deg 

Ax, Ay, A z 

= x, y, and z components of joint translational 
displacement in local body frame, ft 

h 

= Baumgarte error control parameter 

Pa > Pb 

= position vector from point A mass center to 
body A and body B mass center to point B, ft 

co a, co B 

= angular velocity vector of body A and body 
B, relative to inertial reference frame, rad/s 

co x , (o y , ft). 

= body components of inertial angular velocity, 
rad/s 

CO A, CO B 

= angular acceleration vector of body A and 


body B, relative to inertial reference frame, 
rad/s 2 


I. Introduction 

A NALYZING the dynamic separation of multiple bodies within 
the atmosphere is complex and challenging. One problem that 
has received significant attention in the literature is that of store 
separation from aircraft [ J_] . A similar example is the separation of the 
X-15 research vehicle (RV) from the B-52 carrier aircraft [2], In both 
of these cases, the store and the X-15 vehicle are much smaller in size 
than the parent vehicle. The other class of stage separation problem 
involves separation of two vehicles of comparable sizes, as in the case 
of multistage reusable launch vehicles, where the integrity of each 
stage is important after separation. 

NASA studies on stage separation of multistage reusable launch 
vehicles date back to the early 1960s [3-7]. These studies addressed 
the problem of separation of generic two-stage reusable launch 
vehicles. More recently, Naftel and Wilhite [8] and Naftel and Powell 
[9,10] considered staging of dual wing-body vehicles. NASA’s 
interest in stage separation research was renewed in early 2000 in 
connection with design and development of multistage launch 
vehicles. Accordingly, NASA initiated a comprehensive stage 
separation tool development activity that included wind-tunnel 
testing as well as development and validation of computational fluid 
dynamic and engineering-level simulation tools [11] . As a part of this 
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activity, a stage separation analysis and simulation tool called 
ConSep (short for conceptual separation) was developed, which is a 
front end to the commercially available multibody dynamic analysis 
software ADAMS® [ 12 ] . References [13-15] discuss the application 
of ConSep to two-body and three-body separation problems. 
Because it is tied to ADAMS, one disadvantage of ConSep is that it 
cannot be easily integrated into standard trajectory simulation 
software for performing efficient and seamless end-to-end simul- 
ations of launch vehicle trajectories. 

The objective of this paper is to discuss application of the 
constraint force equation (CFE) methodology [ 16 ] for multibody 
stage separation and its implementation in the Program to Optimize 
Simulated Trajectories II (POST2) [T7,^8] for efficient, seamless 
end-to-end simulation of launch vehicle trajectories including stage 
separation. The paper also discusses the verification of the CFE 
methodology and POST2 implementation using three specifically 
designed test cases. Even though this paper discusses the CFE 
implementation in POST2, the method is generic in nature and can be 
implemented in any similar trajectory program capable of simulating 
multiple, unconnected vehicles. 

By itself, POST2 does not have the capability to model internal 
joint forces and moments before separation when the bodies are still 
connected. The CFE implementation in POST2 provides a frame- 
work for computing the internal constraint forces and moments 
acting at joints connecting multiple vehicles. These constraint forces 
and moments can then be applied to each body, together with the 
usual external forces and moments due to gravity, aerodynamics, and 
propulsion. Thus, the CFE methodology simply augments the exter- 
nal loads and does not require modification of the POST2 equations 
of motion. The CFE methodology provides the missing link to model 
accurately the dynamics of generic multibody separation. 

This paper describes the steps taken to verify and validate the CFE 
methodology through its application to three test cases of increasing 
complexity. The verification was performed by comparing results of 
each test case to results obtained from benchmark codes and, when 
possible, flight data. The first test case simulates the motion of two 
rigid bodies connected by a fixed joint, with no external forces. The 
advantage of studying this basic problem is that numerical results can 
be compared with what is to be expected on the basis of first 
principles. The second test case models two rigid bodies connected 
by a universal joint. The third test case consists of a comparison of a 
POST2/CFE simulation of the Mach 7 separation of the X-43A 
vehicle to flight data and to previously published preflight simulation 
results. These three test cases were part of a larger, comprehensive 
series of checks that were performed on the CFE methodology to 
gain confidence in its ability to solve complex constrained motion 
problems using conventional trajectory programs, such as POST2. 
For the first two cases, the CFE simulations are in excellent 


agreement with independent simulations developed using bench- 
mark industry-standard software packages. For the third test case, 
CFE performs as well or better than the preflight simulation, was 
easier to implement, and its results compare reasonably well with the 
flight data. 


II. Constraint Force Equation Methodology 
A. General Formulation 

To illustrate the basic concept of the CFE methodology, consider 
the motion of two rigid bodies connected by a single joint, as shown 
in Fig. L The external forces and moments that act on each body are 
shown in Fig. la. These external forces (F^ EXT) and F^ EXT) ) and 
external torques, or moments of force about mass centers (T® x 1 1 and 
T< EXT) ), are the resultants of gravity, aerodynamic, and propulsive 
forces or moments on each vehicle. Figure fb shows the internal 
constraint forces and moments (F <C0N) and T ,CON) ) acting on each 
vehicle at the joint location. These internal forces and moments 
constrain the way in which one vehicle can move relative to the other, 
and they are dependent upon the external forces acting on each 
vehicle as well as the type of joint. For joints that permit only relative 
rotation between two bodies, as shown in Fig. fb, the forces and 
moments on one body have magnitudes that are equal and directions 
that are opposite to those acting on the other body. The objective of 
the CFE methodology is to compute these internal forces and 
moments and apply them as external forces and moments on each 
body. Figure J_c illustrates the way in which CFE is implemented in 
POST2. At each integration time step, the current version of POST2 
computes the typical external forces and moments acting on each 
vehicle. This information, along with specific geometric information 
about the joint, is provided to the CFE routine, which works in 
parallel to compute the internal forces and moments required to 
satisfy the constraints imposed by the joint. Next, these internal joint 
loads are applied to each vehicle as additional external forces and 
moments, and the POST2 solution is propagated in the usual manner 
to the next time step. Thus, the net external forces and moments on 
each vehicle are the sum of 1 ) the usual external forces and moments 
and 2) the joint loads applied to each vehicle as additional external 
forces and moments. Consequently, the CFE joint model simply 
augments the vehicle external loads and does not require 
modification to the POST2 equations of motion. 

The equations of constrained motion of two rigid bodies (A and B ) 
connected by a single joint (point A in body A and point B in body B) 
are as follows. For body A, 

Ff^ + Pj CON) = m A x A (1) 




a) External forces and moments b) Internal forces and moments c) Resultant forces and moments 
Fig. 1 Illustration of basic concept of CFE methodology. 
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Tf^ + p A x F< con > + T<“ N) = + (2) 

where p A is the position vector from the mass center of A to point A of 
A at which the constraint force is applied. Similarly, for body B, 

F ' RXT) + F r N ) = mRXR (3) 

Tr T) + Pb x F< CON) + Ti C0N) = l B -io B + a> B xl B -co B (4) 

Note that, in the case of unconstrained (free) motion, the constraint 
forces and moments are equal to zero, and Eqs. (j_ — 4) reduce to the 
standard equations of motion for two free rigid bodies. Each mass 
center acceleration, rigid body angular acceleration, constraint force, 
and constraint torque are regarded as unknowns; therefore, Eqs. (1— 
4) represent 12 scalar equations in 24 unknowns. Hence, 12 
additional scalar equations are required. Of these, six are obtained 
from two vector relationships that are consequences of the law of 
action and reaction: 

F (CON) + F (CON) = 0 (5) 

T f 0N) + Xr N) + (r« - r A ) x Fg C0N) = 0 (6) 

where r A is the position vector from a point fixed in an inertial 
reference frame to A. Because x A is the position vector from a point 
fixed in an inertial reference frame to the mass center of A, we can 
write r A = x A + p A . Note that vector p A is fixed in A. Position 
vector is defined similarly. 

The six remaining equations fall into one of two categories. First, 
equations of constraint are formed that describe how a particular joint 
restricts relative translation or relative rotation. Second, by 
considering the translation or rotation that is permitted by the joint, 
equations are written to account for the absence of constraint force or 
torque in a particular direction in view of the ideal (perfectly smooth) 
nature of the joint. 

When relative translation is constrained, the distance between two 
points must remain fixed in a particular direction. That is, 

(r B - r A ) • e A = 0 (7) 

where the joint prevents relative translation in the direction of unit 
vector e A , which is fixed in body A. One equation having the form of 
Eq. (7) is needed to account for each direction in which translation is 
constrained. For example, a fixed joint constrains translation in three 
orthogonal directions; therefore, three equations in the form of 
Eq. (7) are required. In each such equation, the role of e A is played by 
one of three mutually orthogonal unit vectors fixed in A. For a 
prismatic (sliding) joint that permits translation in only one direction 
(and restricts translation in two perpendicular directions), two 
equations having the form of Eq. (7) are required. 

A constraint on relative rotation can be viewed as a requirement 
that two unit vectors must remain perpendicular; each unit vector is 
fixed in one of the bodies and is normal to the axis about which 
rotation would take place if the constraint were not present. The 
constraint is expressed by setting the scalar product of the two unit 
vectors equal to zero: 

e rj ' e A = 0 (8) 

where the two unit vectors, e A and e fi (one fixed in each body), are 
chosen to be perpendicular to each other throughout the constrained 
motion. One equation in the form of Eq. (8) is required for each 
direction about which rotation is constrained. 

In their current form, Eqs. (7) and (8) are not coupled with Eqs. (J- 
4). For this purpose, the constraint equations (7) and (8) are 
differentiated twice with respect to time in an inertial reference frame 
so that the resulting equations involve the unknown linear and 
angular accelerations of bodies A and B, and thus can be coupled to 


the equations of motion. To this end, the second derivatives that need 
to be evaluated are 

-^2 — f a) ' e A = 0] (9) 

d t 2 

^[e B -e A = 0] (10) 

df 2 

The result of differentiating Eq. (7) once is 

(r B ~ r A ) • e A + (r B - r A ) • (« A x e A ) = 0 (11) 

where r A and r B are, respectively, the time derivatives in an inertial 
reference frame of r A and r B . Differentiating once more, one obtains 

(x s + x p B - x A - u> A x p A ) ■ e A + (r B - r A ) • (<a A x e A ) 

= 2(r fl - r A ) • (e A x co A ) - (r B - r A ) • (t» A x u> A x e A ) 

+ [("a x a A x p A ) - (a> B x (o B x p B )] • e A (12) 

Differentiating once with respect to time in Eq. (8) yields 

{a> B x e B ) • e A + e fl • (« A x e A ) = 0 (13) 

and differentiating once more gives 

(d) B - d) A ) • (e B x e A ) = (h) B - fc> A ) • [e A x (u> B x e B ) 

- e B x (<u A x e A )] (14) 

When relative translation is permitted in a certain direction, an 
equation having the form of Eq. (_12) is not required. Likewise, when 
relative rotation is allowed, a relationship in the form of Eq. (_14) is 
not necessary. Instead, the required equations can be obtained using 
the condition that the constraint forces or moments are zero in those 
directions (with the assumption that the surfaces of the joint are 
perfectly smooth). Thus, for a joint that permits translation in a 
certain direction, 

F (CON). e = 0 (15) 

Depending on the circumstances, F <C0N) can be the constraint force 
acting either on A or on B. Similarly, when a joint permits relative 
rotation about a certain direction, 

T (con) . e _ o ( 16 ) 

In summary, equations having the form of Eqs. (1_2) and (14-16) 
provide, in combination, a total of six scalar equations. 

B. Joint Stabilization 

Equations ( 1-16) represent a sufficient set for solving the general 
constraint load problem described herein. However, the accuracy of 
the constraint load solution is sensitive to computational numerical 
error and initial joint misalignment. In particular, an accumulation of 
numerical error typically manifests itself as joint separation and/or 
misalignment error over time. This in turn introduces additional 
errors in the solution, and the cycle continues, leading to a buildup of 
error over time. To reduce, or even reverse, joint separation and 
misalignment errors, the CFE algorithm implements an optional 
stabilization technique known as Baumgarte [J_9] stabilization. 
Ordinarily, the constraint equations ( 12 ) and (14) take the following 
generalized form: 

g = 0 (17) 

where g represents either of the (nondifferentiated) constraints in 
Eqs. (7) and (8). The stabilization technique is implemented by 
augmenting Eq. (17) with terms involving the once differentiated and 
nondifferentiated forms of g: 

g + 2t?g + Tg = 0 (18) 
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Thus, the new constraint relation in Eq. (J_8) considers not only the 
acceleration constraints but also the corresponding velocity and 
original constraints. These additional terms serve as a form of 
proportional and derivative control that can mitigate joint separation 
and misalignment errors. Since the stabilization effect is dependent 
on the specific dynamic system, the Baumgarte parameter )/ can be 
tuned to control the degree of constraint error for a particular 
application. 



Fig. 2 Testcase 1: two rigid bodies.A and /i. connected by afixed joint. 


C. Implementation 

To implement the CFE method for a single joint connecting two 
rigid bodies, 24 scalar equations are formed from Eqs. (l_-6), (12), 
and (14-16) that are linear in 24 unknowns. These equations can be 
expressed in matrix form, Ax = b, and solved accordingly using a 
standard matrix inversion technique. The column matrix x contains 
the 24 unknown parameters, which are the scalar components of the 
following vectors: x A , x B , ch A , <i> B , F A CON) , F B CON) , T A CON) , and 
T< CON ’. Matrix A can be singular when, for example, one neglects a 
central principal moment of inertia, such as for a slender rod, or when 
the mass of a body is neglected. The matrix also becomes singular 
when constraint equations are redundant; in other words, when one is 
linearly dependent on the others. An example of a redundant 
constraint equation would be a body connected by two joints that 
both constrain translation in a given direction. Thus, care should be 
taken to formulate the problem in a way that avoids these difficulties. 
The X-43A test case discussed in Sec. 1 II.C explains how the 
problem’s formulation was altered in order to avoid redundant 
constraint equations. 

In general, the CFE method can be applied to problems containing 
multiple bodies and multiple joints, so long as there are not redundant 
constraints. In general, the number N of equations and unknowns 
involved in a problem is given by 

N = 6 (N b + Nj) (19) 

where N B is the number of rigid bodies, and 2 Nj is the number of 
joints. 

Finally, it is possible to improve computational efficiency by 
implementing the CFE algorithm in a way that reduces the number of 
equations and unknowns. One technique that has been successfully 
employed by the authors is the conventional minimal matrix 
approach [20-22], This technique combines Eqs. O), (2), (12), and 
(14) to solve only for the unknown constraint loads rather than 
solving for N unknowns. The minimal matrix approach is 
advantageous, since it involves only those constraint equations that 
are needed for any particular joint configuration; thus, no extraneous 
calculations are made. 

III. Test Cases 

A suite of test cases was developed to check and verify the CFE 
methodology for application to multibody separation problems in the 
launch vehicle staging environment. These test cases are listed in 
Table J_ and cover the types of joints typically needed to model stage 
separation problems. The list is not exhaustive and, in general, it is 
possible to model joints that simultaneously constrain any 


Table 1 Test cases performed for CFE verification 
(2-D denotes two-dimensional) 


Case 

Joint type 

Translational 
degrees of 
freedom 

Rotational 
degrees of 
freedom 

i 

Fixed 

0 

0 

2 

Revolute (hinge) 

0 

1 

3 

Universal 

0 

2 

4 

Ball 

0 

3 

5 

1 -D translational (slider) 

1 

0 

6 

2-D translational (planar) 

2 

0 

7 

Cylindrical 

1 

1 


combination of translational and rotational degrees of freedom with 
CFE simply by forming the appropriate equations described in 
Eqs. (12) and (14-16). 

A description of the use of CFE to solve test cases 2-7 from Table 1_ 
can be found in [23], Two test cases from Table J_ (the fixed joint and 
the universal joint) are discussed in this paper to illustrate the testing 
and verification that was performed on the CFE methodology. A third 
test case was conducted that used POST2/CFE to model the X-43A 
stage separation problem as a one-dimensional (1-D) translational 
(or sliding) joint in series with a revolute (or hinge) joint. 

A. Test Case 1: Fixed Joint 

This test case involves two rigid bodies, denoted as bodies A and 

B, connected by a fixed joint (see Fig. 2). The mass properties for 
each body are shown in Table 2. The fixed joint constrains point A 
attached in A to remain coincident with point B attached in B. Three 
constraint equations having the form of Eq. (7) can be written as 

( r B ~ r A ) • a r = 0 (r = 1,2,3) (20) 

where the role of e A inEq. (7) is played in turn by a [, a 2 , and a 3 , each 
of which belongs to a set of three right-handed mutually orthogonal 
unit vectors fixed in A, as shown in Fig. 2. Suppose that b,, b 2 , andb 3 
are a similar set of unit vectors fixed in B, such that b r has the same 
direction as a r (r = 1,2,3) when A and B are attached to each other. 
In that case, the fixed joint constrains certain unit vectors fixed in A to 
remain perpendicular to other unit vectors fixed in B. For example, 
one can write the following three constraint equations having the 
form of Eq. (8): 


= 0 

(21) 

= 0 

(22) 

= 0 

(23) 


Thus, a total of six constraint equations have been written. As 
mentioned previously, these constraint equations are differentiated 
twice with respect to time so that the resulting relationships involve 
the unknown accelerations and angular accelerations. Thus, one 
employs three equations having the form of Eq. (12), in which the 
role of e A is played by a,, a 2 , and a 3 , respectively. Likewise, three 
equations having the form of Eq. (J_4) come into play. In the first of 
these, the roles of e A and e B are played by a ( and b 2 . In the second 
such equation, a t and b 3 play the parts of e A and e B . In the third 
relationship, a 3 and b 2 are substituted for e A and e B . 

No relationships having the form of Eq. (15) or Eq. ( 16 ) are 
applicable in the case of the fixed joint, because no relative motion is 
permitted by the joint. 


Table 2 Mass properties for test case 1, fixed joint 



Body A 

Body B 

m 

1 slug 

4/10 slug 

I XX 

5/12 slug ■ ft 2 

1 /6 slug • ft 2 


26/12 slug • ft 2 

1 /6 slug ■ ft 2 

hz 

29/2 slug • ft 2 

4/15 slug • ft 2 
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Fig. 3 Test Case 2: Two rigid bodies, A and B, connected by a universal 
joint. 

B. Test Case 2: Universal Joint 

This test case involves two rigid bodies, A and B , for which the 
mass is distributed uniformly, and they are connected by a universal 
joint. The vector from the joint to A* is— 4.854a! — 0.096a 2 , and the 
vector from the joint to B* is 3.851b, — 0.096b 2 . This configuration 
is shown in Fig. 3. The universal joint constrains point B fixed in B to 
remain coincident with point A fixed in A at the intersection of the 
two arms in the cross. Therefore, the three translational constraint 
equations presented in Eq. (20) for the fixed joint apply in this case as 
well. 

For this problem, one rotational constraint equation having the 
form of Eq. (8) is needed. If a 2 is chosen to be parallel to the ai m that 
is fixed in A, and b 3 is chosen to be parallel to the arm that is fixed in 

B , then the rotational constraint equation imposed by the universal 
joint can be described with the relationship 

b 3 • a 2 = 0 (24) 

The foregoing four constraint equations are differentiated twice with 
respect to time. Thus, for the translational constraint, one employs 
three equations having the form of Eq. 02), in which the role of e A is 
played by a,, a 2 , and a 3 , respectively. Likewise, for the rotational 
constraint, the result is obtained by employing Eq. (14) and making 
two substitutions: a 2 for e A , and b 3 for e fl . 

The other two equations are obtained with the aid of Eq. (16): it is 
assumed that the universal joint is ideal (frictionless), so that no 
constraint torque is exerted in the directions of the arms in the cross: 

T A C0N) • a 2 = 0. T<! CON) • b 3 = 0 (25) 

The mass properties of each body are summarized in Table 3 
(products of inertia are zero for each body in this problem). 

C. Test Case 3: X-43A 

This test case was selected to evaluate the capability of POST2/ 
CFE to model and simulate a realistic stage separation problem by 
comparing results with flight data from NASA's X-43A hypersonic 
scramjet test vehicle (see Fig. 4). In 2004, NASA conducted two 
successful scramjet test flights of the X-43A RV at speeds near 
Mach 7 and 10. The roughly 10 s scramjet test was performed after 
the RV was boosted to the target flight condition by the Hyper-X 
launch vehicle (HXLV). A dimensioned drawing of the mated RV 
and HXLV is shown in Fig. 5. Additional information describing the 
X-43A RV and HXLV may be found in [24]. 

The X-43 A stage separation event was initiated approximately 3 s 
after HXLV burnout at a dynamic pressure of 1000 psf. Separation 
began when two pyrotechnically actuated pistons extending from the 
booster pushed against the RV to induce ~ 1 6 ft/s of relative velocity 
between the two vehicles. The two pistons, which were positioned 


Table 3 Mass properties for test ease 2, 
universal joint 



Body A 

Body B 

m 

622.0 slugs 

622.0 slugs 

I XX 

220.3 slug ■ ft 2 

881.2 slug -ft 2 


3524.7 slug ■ ft 2 

2695.3 slug ■ ft 2 

! zz 

3330.3 slug ■ ft 2 

1917.8 slug -ft 2 



Fig. 4 Test case 3: artistic rendering of X-43A (RV) separation from 
HXLV booster. 



Fig. 5 Test case 3: dimensioned drawing of mated X-43A RV and 
HXLV. 


roughly 9 in. on each side of the RV centerline, were initially in 
contact with a cuplike ball joint attached to the RV that permitted 
rotation about the piston contact point and translation only along the 
line of action of the piston force, which was oriented 4 deg below the 
horizontal. Both pistons were connected to the same gas chamber and 
were activated at the same time and remained in contact with the RV 
for ~0. 1 sec until they reached the end of their 9 in. stroke length. 
Although the ball joint permitted rotation in any direction at each 
individual piston contact point, the fact that both pistons were fired 
concurrently and were simultaneously in contact with the RV 
effectively excluded relative rolling and yawing rotation between the 
two separating vehicles. 

Once the pistons reached the end of their 9 in. extensions, the RV 
separated from the HXLV. At this time, the HXLV commanded a 
nosedown pitching moment to move away from the RV, and the RV 
was steered to its required target condition. 

To test the CFE methodology, this problem was set up using 
POST2/CFE to model the piston constraint forces. For the purpose of 
this paper, the two pistons were modeled as a single piston to avoid 
overconstraining the problem. Thus, a single piston was modeled as a 
sliding joint in series with a revolute joint; that is, relative translation 
was allowed in one direction (the piston line of action) and relative 
rotation was permitted only about the pitch axis. The piston was 
placed at the centerline of the RV, and the axial force that pushed the 
vehicles apart was doubled. A schematic highlighting these 
modeling details is shown in Fig. 6. 

To illustrate how the constraint equations for this joint are derived, 
unit vectors are defined in the HXLV such that a, is the direction in 
which translation or sliding is permitted (along the piston line of 
action) and a 2 indicates the direction in which rotation is free to occur 
(pitch axis). A similar set of unit vectors b r can be defined that are 
fixed in the RV and initially aligned with the unit vectors a r fixed in 
the HXLV (see Fig. 6). For this case, two constraint equations having 
the form of Eq. (7 ) can be written as 

(ffl — r A ) • a r = 0 (r = 2,3) (26) 

As before, these two constraint equations are differentiated twice 
with respect to time, as in Eq. (12). Now, a 2 and a 3 are substituted for 
e A . Furthermore, two constraint equations describing the restriction 
imposed on relative orientation are formed: 
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Fig. 6 Test case 3: CFE modeling of piston contact for X-43A 
separation problem. 

b i • a 2 = 0. b 3 • a 2 = 0 (27) 

Similarly, these constraint equations are differentiated twice with 
respect to time, as in Eq. ( _14), where a 2 is now substituted for e A , and 
b, and b 3 are substituted for e B . With the assumption that the joint is 
frictionless, the constraint force and moment in the directions where 
motion is permitted become, using Eqs. (15) and (16), 

F <j CON) • a, = 0, T<j CON) • a 2 = 0 (28) 

IV. Results and Discussion 

A. Test Case I: Fixed Joint 

The two rigid bodies, A and B, are assumed to be rigidly connected 
to each other for 10 s (see Fig. 2) and then released instantaneously. 
No external forces or moments were assumed to act on either body 
for the entire duration of the 20 s simulation. Initial conditions 
correspond to zero inertial velocity of the system mass center, S* , and 
an inertial angular velocity of the rigid body formed by A and B 
equal to 

63.02aj — 82.32a 2 + 80.25a 3 deg /s 

After 10 s, the joint connecting A and B is released, and their motions 
become unconstrained. Because there are no external forces or 
moments, integrals of the motion involving linear and angular 
momentum must remain constant, and this fact can be used to check 
numerical solutions of the equations of motion. Therefore, one can 
expect the simulation results to show that, after separation, the mass 
centers of A and B travel in straight lines with constant velocities. 

To verify the POST2/CFE simulation results, a simulation was 
also created with AUTOLEV [25-27], an interactive program 
designed specifically for the kinematic and dynamic analysis of 
mechanical systems. Symbolic manipulation (computer algebra) is 
used to formulate explicit equations of motion for the particular 
system of interest using any method; the software is particularly well 
suited to Kane’s method [27,28], For this test case, AUTOLEV was 
used to create a computer program to simulate the motions of A and B 
and determine the constraint forces required to hold them together for 
the first 10 s. 

The AUTOLEV results were generated using a variable-step 
integrator with an absolute error limit of 1 x 1CL 8 and a relative error 
limit of 1 x 10~ 7 . The POST2/CFE method employed a fixed-step 
integrator with a step size of 0.0001 s. The CFE routine also used a 
Baumgarte parameter of 5.0. 

As expected, both linear and angular momentums remained 
constant throughout the simulation. Time histories for the inertial x, 
y, and z components of velocity of the mass center of body A are 
shown for 20 s in Fig. 7. The values oscillate during the first 10 s, after 
which the joint is released and their values become constant. Figure 7 
shows that there is an excellent agreement between POST2/CFE and 
AUTOLEV results. Similarly, the velocity components for the mass 
center of body B in Fig. 8 also indicate an excellent agreement 
between POST2/CFE and AUTOLEV. 



Fig. 7 Test case 1, body A: velocity components of mass center. Joint 
release occurs at t = 10 s (circles: POST2/CFE; lines: AUTOLEV). 



Fig. 8 Test case 1, body B : velocity components of mass center. Joint 
release occurs at t = 10 s (circles: POST2/CFE; lines: AUTOLEV). 

The angular velocity time history comparisons for bodies A and B 
are shown in Figs. 9 and H), respectively. Again, there is excellent 
agreement between the two approaches. The change in the angular 
velocity of body B at the time of release is more noticeable than for 
body A, because B is smaller and has a mass distribution that differs 
from that of the composite body before joint release. The main results 
of the CFE method are the constraint forces and moments, which are 
shown in Figs. _H and 12 and are applied as additional external loads 
within POST2. The figures show that the forces and torques 
determined by CFE agree very well with the AUTOLEV results. 



Fig. 9 Test case 1, body A: angular velocity components. Joint release 
occurs at t = 10 s (circles: POST2/CFE; lines: AUTOLEV). 
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Fig. 10 Test case 1, body B: angular velocity components. Joint release 
occurs at t = 10 s (circles: POST2/CFE; lines: AUTOLEV). 



Fig. 13 Test case 1: joint displacement in x, y, and z directions for 
POST2/CFE method during time when bodies are connected. Joint 
release occurs at t = 10 s. 



Fig. 11 Test case 1: comparison of constraint forces required to hold 
bodies together. Joint release occurs at t = 10 s (circles: POST2/CFE; 
lines: AUTOLEV). 


An important metric to assess the POST2/CFE methodology is the 
relative joint displacement between the two bodies when they are 
supposed to stay connected. This parameter is computed as the 
position vector from A to B and should be zero while the joint 
constraint is imposed for the first 10 s. Hence, any deviation from 
zero is a measure of accuracy of the CFE algorithm. In the CFE 
algorithm, the constraints are satisfied identically at the acceleration 
level [Eqs. (12) and (14)]; however, as with algorithms of this type. 



Fig. 12 Test case 1: comparison of constraint torques required to hold 
bodies together. Joint release occurs at t = 10 s (circles: POST2/CFE; 
lines: AUTOLEV). 


the constraints at the velocity and position levels are inevitably 
subject to numerical integration errors. Such errors are functions of 
the step size in a fixed-step integration scheme or the error limits in a 
variable-step approach. Baumgarte stabilization controls these errors 
from growing arbitrarily large during a simulation. 

The choice of step size and Baumgarte factor r) for this problem 
represents a reasonable balance between CPU time and error buildup. 
Figure L3 shows the joint translational displacement as a function of 
time for the 10 s period when the bodies are supposed to stay 
connected. In this paper, the quantities Ax, Ay, and Az are the x, y, 
and z components of joint translational displacement in the local 
body frame, and they are defined as follows: 

Ax = (I# - r A ) • a, , Ay = (r„ - r A ) • a 2 

Az = (r B - r A ) • a 3 (29) 

In the interest of clarity, only the POST2/CFE results are shown. The 
magnitude of this displacement distance is below 0.01% of the total 
body length. 

B. Test Case 2: Universal Joint 

To check and verify POST2/CFE results, an independent 
simulation using ADAMS software was generated. At t = 0, each 
body was assumed to be at rest. A time-varying force is applied to 
body A at a location relative to the center of mass given by the vector 
— 2.0a£ + 0.5a 3 ft. The components of this time-varying force are 
shown in Fig. f4. The mass properties of each body are summarized 
in Table 3. The products of inertia are zero for each body in this 
problem. 



Fig. 14 Test case 2: external force components applied to body A 

^p(EXT) _ p(EXT)^ 
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The motion of bodies A and B was simulated for a duration of 5 s. 
The POST2/CFE used a fixed-step-size integrator with a step size of 
0.0001 s and a Baumgarte factor, r] = 5.0. The results of the 
simulations are shown in Figs. 15-18. The angular velocity of body A 
is shown in Fig. 15. The constraint force and moment components are 
shown in Figs. 1T> and 17, respectively. There is excellent agreement 
between POST2/CFE and the ADAMS solutions, as evidenced by 
these figures. Finally, the relative joint translational displacement 
time histories are shown in Fig. j_8. As in the previous test case of the 



Fig. 15 Test case 2: body B: angular velocity components (circles: 
POST2/CFE; lines: ADAMS). 



Fig. 16 Test case 2, body B : Comparison of constraint force (circles: 
POST2/CFE; lines: ADAMS). 



Fig. 17 Test case 2, body B: Comparison of constraint torque (circles: 
POST2/CFE; lines: ADAMS). 


fixed joint problem, the joint displacement is relatively small, well 
below 0.005% of body length. 

C. Test Case 3: X-43A 

The objective of the third test case was to evaluate the capability of 
POST2/CFE to model a realistic stage separation problem. For the 
X-43 A problem that was chosen, POST2/CFE results of the angle of 
attack and angle of sideslip during separation were compared with 
data recorded during flight [29], In addition, because many of the 
parameters needed to fully verify the POST2/CFE results were not 
included in the flight data, further comparisons were made with 
previously published simulation results that were part of the X-43A 
preflight analysis [30] . This preflight simulation, which was also 
developed using POST2, modeled each joint between the two pistons 
and the RV as a massless spring, with a prescribed force profile in the 
direction along the piston centerline, and with restoring forces 
computed from beam deflection theory in the lateral directions 
perpendicular to the axis of piston motion. By modeling the pistons 
as massless springs, the need to consider constraint equations in the 
POST2 formulation was avoided; however, the spring-model 
implementation introduced modeling complexity and required 
numerous customizations to the POST2 source code. 

Reference [ 30 ] details the development of this preflight POST2 
spring-model simulation and describes the excellent level of 
agreement between its results and the results of an even higher 
fidelity simulation of the X-43A stage separation that was developed 
using ADAMS. Because of this heritage, the POST2 spring-model 
simulation from [ 30 ] was considered suitable for verifying the 
POST2/CFE simulation results. 

A description of how the X-43 A stage separation problem was set 
up in the POST2/CFE formulation was presented in Sec. III. Because 
of the generalized implementation of the CFE method within 
POST2, the problem was much easier to set up than the spring-model 
implementation, and it required few modifications to the POST2 
source code. The primary difference in problem setup between the 
two approaches was in modeling of the pistons, where POST2/CFE 
modeled only one piston to avoid overconstraining the problem 
while the spring-model implementation modeled both pistons. Since 
the objective of the simulation was to determine the motion of both 
vehicles during separation to assess vehicle recontact potential and 
optimize control system settings, the details of how the pistons were 
modeled were not important as long as the motion was constrained 
the same way while both vehicles were joined. As a measure of how 
well the joint constraint was satisfied, the relative joint displacement 
between the two separating vehicles is shown in Fig. 19. This 
parameter is the distance from the RV joint location to the F1XLV 
joint location, and it was computed in the piston axial (line of action, 
a ( ), lateral (a 2 ), and normal (a 3 ) directions. The results indicate that 
CFE does a better job of enforcing the translational constraint, and 
there is virtually no joint displacement in the lateral or normal 



Fig. 18 Test case 2: joint displacement in x, y, and z directions for 
POST2/CFE method during time when bodies are connected. 
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Fig. 19 Test case 3: comparison of joint displacement in x, y, and z 
directions for CFE and spring model. 
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Fig. 20 Test case 3: comparison of relative angular displacement 
between POST2/CFE and spring model. 


directions. Similar results are seen in Fig. 20, which shows the 
relative difference in attitude between the two vehicles. Again, CFE 
does an excellent job of enforcing the rotational constraint, 
essentially eliminating any relative roll and yaw between the RV and 
HXLV. The constraint forces and moments about the vehicle center 
of mass that were computed by each technique are shown in Figs. 2J_ 
and 22, respectively. In these plots, the forces and moments 



Time, s 

Fig. 21 Test case 3: comparison of constraint forces computed by 
POST2/CFE and spring model. 
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Fig. 22 Test case 3: comparison of constraint moments computed by 
POST2/CFE and spring model. 


computed by the spring model include the contributions of both 
pistons. The results show that both the spring model and POST2/CFE 
compute very similar forces and moments. Flowever, the spring- 
model results are oscillatory in nature, since the joint constraint is not 
satisfied exactly. 

Flight data comparisons of the angle-of-attack and angle-of- 
sideslip profiles are shown in Fig. 23. Both the spring-model and 
POST2/CFE results match very well and are nearly identical in the 
angle-of-attack response. However, there are slightly larger differ- 
ences in the sideslip response. In addition, both simulation 
techniques match the observed flight data well. It should be noted that 
the degree to which the simulation results match the flight data 
depends on other factors besides CFE modeling, such as aerody- 
namic force and moments, winds, accuracy of flight data instrumen- 
tation, etc. Aerodynamic modeling was especially critical because 
separation occurred at a dynamic pressure of 1000 psf. Thus, this 
problem is much more complex than the first two test problems where 
the external forces and moments were prescribed. 

To illustrate the sensitivity of the simulation results to aerody- 
namic modeling, computations were made using the nominal 
(preflight) aerodynamic database and a postflight reconstructed 
aerodynamic database in which the uncertainty terms had been 
adjusted to match the flight data more closely. The variations in angle 
of attack and sideslip differ significantly depending on which 
aerodynamic model was used. Note that all POST2/CFE results were 
generated using the adjusted aerodynamics. 

As a further comparison, results are also shown for an 
unconstrained case where no constraint forces were applied to model 
the joint (but the piston axial force that pushed the vehicles apart was 
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Fig. 23 Test ease 3: simulation vs flight data comparison of angle-of- 
attack and angle-of-sideslip profiles (aero: aerodynamic, adj: adjusted). 
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Fig. 24 Test case 3: simulation vs flight data comparison of x, y, and z 
accelerations in local body frame. 

still modeled). During the first second of separation, the 
unconstrained results differed significantly from the cases that 
included the constraint, demonstrating the importance of accurate 
modeling of the piston joint. Even though the differences between the 
constrained and unconstrained angle-of-attack and sideslip profiles 
disappeared after the control system became active and was able to 
compensate, it is important to capture the dynamics during the first 
second of flight, since that was when the risk of recontact was 
highest. 

Similar comparisons are shown for the linear accelerations in 
Fig. 24. The vertical (i.e., body z direction) acceleration results were 
the most sensitive to the way in which the piston constraint was 
modeled, and they were the most important for assessing the risk of 
recontact. Again, the spring-model and CFE results were nearly 
identical and provided a good match to the flight data. The 
differences between the constrained and unconstrained cases were 
actually larger than the differences due to aerodynamics. 

Figure 25 compares the angular velocity results with the flight 
data. Once again, there is excellent agreement between the spring- 
model and CFE simulation results. Both simulations match the pitch 
rate flight data very well and do a much better job at predicting the 
pitching motion than the unconstrained case. The comparison 
between simulation and flight for yaw rate and roll rate is not as good. 
Much of the difference can be attributed to inadequate lateral- 
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Fig. 25 Test case 3: simulation vs flight data comparison of body roll, 
pitch, and yaw rate profiles. 


directional aerodynamic modeling, as is evident in the large 
difference in results between the cases with nominal aerodynamics 
and adjusted aerodynamics. Again, capturing the pitch plane 
dynamics was important for assessing the risk of recontact. 

V. Conclusions 

The CFE methodology provides a simple method to calculate the 
internal constraint forces and moments acting at joints connecting 
multiple bodies, and it applies them as external forces and moments. 
The details of the CFE methodology and its implementation in the 
POST2 are discussed in this paper. The CFE methodology is generic 
in nature and can be implemented in any trajectory simulation 
software. To illustrate the capability of the CFE algorithm, three test 
cases of increasing complexity were presented. The first case 
examined the motion of two force-free rigid bodies connected by a 
fixed joint, subject to an initial rotation. The simple nature of this 
problem made it possible to compare simulation results to analytical 
predictions based on first principles. In addition, simulating a fixed 
joint presented a challenging case with which to assess the relative 
displacement between the two bodies at the joint location, a key 
performance metric for the CFE algorithm. The second case involved 
rigid bodies connected by a universal joint, and it typified the testing 
that was performed on a number of different joints applicable to stage 
separation problems. The CFE results were in excellent agreement 
with results obtained using AUTOLEV in the first test case and with 
those obtained using ADAMS in the second case. 

In the third test case, the ability of the CFE method to model 
realistic launch vehicle stage separation was demonstrated by 
simulating X43-A stage separation. These results agreed well with 
previously published results that were based on an engineering 
model, and they actually performed better at satisfying the joint 
constraint. Moreover, the CFE results agreed reasonably well with 
the flight data, and it was shown that the level of agreement degraded 
significantly if the constraint was not modeled. 

The success of these test cases provides confidence in the ability of 
the CFE methodology to solve constrained motion problems in 
conventional trajectory simulation programs and to provide the 
capability of performing generic end-to-end simulations of launch 
vehicle trajectories, including stage separation. 

Acknowledgments 

The authors gratefully acknowledge Anne Rhodes for her 
invaluable assistance in preparing this paper for publication and Jim 
Beaty for his helpful comments. 

References 

[1] Dillenius, M. F. E., Perkins, S. C., and Nixon, D., "Pylon Carriage and 
Separation of Stores ,” Tactical Missile Aerodynamics-General Topics, 
edited by M. J. Hemsch, Progress in Astronautics and Aeronautics, 
AIAA, Washington, D.C., 1992, pp. 575-666. 

[2] Taylor, R. T., and Alford, W. J., Jr., "A Wind Tunnel Investigation of the 
Carry Loads and Mutual Interference Effects of 1/40-Scale Models of 
the X-15 and B-52 Airplanes In Combination,” NASA TM X-184, 
1959. 

[3] Decker, J. P, and Wilhite, A. W., “Technology and Methodology of 
Separating Two Similar Size Aerospace Vehicles Within the 
Atmosphere,” AIAA Paper 1975-0029, Jan. 1975. 

[4] Decker, J. P., “Experimental Aerodynamics and Analysis of the Stage 
Separation of Reusable Launch Vehicles,” NASA, SP-148, 1967. 

[5] Decker, J. P., and Gera, J.. "An Exploratory Study of Parallel-Stage 
Separation of Reusable Launch Vehicles,” NASA TN D-4765, 1968. 

[6] Decker, J. P., "Aerodynamic Interference Effects Caused by Parallel- 
Staged Simple Aerodynamic Configuration at Mach Numbers of 3 and 
6,” NASA TN D-5379, 1969. 

[7] Wilhite, A. W., "Analysis of Separation of the Space Shuttle Orbiter 
from a Large Transport Airplane,” NASATM X-3492, 1977. 

[8] Naftel, J. C., and Wilhite, A. W., “Analysis of Separation of a Two-Stage 
Winged Launch Vehicle,” AIAA Paper 1986-0195, Jan. 1986. 

[9] Naftel, J. C., and Powell, R. W., “Aerodynamic Separation and 
Glideback of a Mach 3 Staged Orbiter," AIAA Paper 1990-0223, 
Jan. 1990. 


TARTABINI ET AL. 


583 


[10] Naftel, J. C., and Powell, R. W., “Analysis of the Staging Maneuver and 
Booster Glideback Guidance for a Two-Staged, Winged, Fully 
Reusable Launch Vehicle,” NASA TP-3335, 1993. 

[11] Murphy, K. J., Buning, P. G., Pamadi, B. N., Scallion, W. I., and Jones, 
K. M., “Status of Stage Separation Tool Development for Next 
Generation Launch Vehicle Technologies,” ALAA Paper 2004-2595, 
2004. 

[12] Using ADAMS/Solver, Mechanical Dynamics, Ann Arbor, MI, 1999. 

[13] Pamadi, B. N., Neirynck, T. A., Coveil, P. F., Hotchko, N. J., and Bose, 
D. M., “Simulation and Analyses of Staging Maneuvers of Next 
Generation Reusable Launch Vehicles,” ALAA Paper 2004-5 1 85, 2004. 

[14] Pamadi, B. N., Neirynck, T. A., Hotchko, N. J., Scallion, W. I., Murphy, 
K. J., and Covell, P. F., “Simulation and Analyses of Stage Separation of 
Two-Stage Reusable Launch Vehicles,” Journal of Spacecraft and 
Rockets , Vol. 44, No. 1, 2007, pp. 66-80. 

doi: 10.2514/l. 17896 

[15] Pamadi, B. N., Hotchko, N. J., Samareh, J., Covell, P. F., and Tartabini, 
P. V., “Simulation and Analyses of Multi-Body Separation in Launch 
Vehicle Staging Environment,” AIAA Paper 2006-8033, 2006. 

[16] Toniolo, M. D., Tartabini, P. V., Pamadi, B. N., and Hotchko, N. J., 
“Constraint Force Equation Methodology for Modeling Multi-Body 
Stage Separation Dynamics,” AIAA Paper 2008-0219, 2008. 

[17] Brauer, G. L., Cornick, D. E., and Stevenson, R., “Capabilities and 
Applications of the Program to Optimize Simulated Trajectories 
(POST),” NASA CR-2770, 1977. 

[18] Striepe, S. A., Powell, R. W., Desai, P. N., Queen, E. M., Brauer, G. L., 
Cornick, D. E., Olson, D. W., Petersen, F. M., Stevenson, R., Engel, M. 
C., Marsh, S. M., and Gromko, A. M., Program to Optimize Simulated 
Trajectories ( POST II), Vol. 2: Utilization Manual, Ver. 1.1.7, NASA 
Langley Research Center, Hampton, VA, July 2002. 

[19] Baumgarte, J., “Stabilization of Constraints and Integrals of Motion in 
Dynamical Systems,” Computer Methods in Applied Mechanics and 
Engineering, Vol. 1, No. 1, 1972, pp. 1—16. 

doi: 1 0. 1 0 1 6/0045-7825(72)900 1 8-7 

[20] Greenwood, D. T., Principles of Dynamics, Prentice-Hall, Englewood 


Cliffs, NJ, 1965, pp. 267-269. 

[21] Wang, J. T., and Huston, R. L., “Kane’s Equations with Undetermined 
Multipliers: Application to Constrained Multibody Systems,” Journal 
of Applied Mechanics, Vol. 54, No. 2, 1987, pp. 424^4-29. 

doi: 10.1 1 15/1.3173031 

[22] Anderson, K. S., “An Order n Formulation for the Motion Simulation of 
General Multi-Rigid-Body Constrained Systems,” Computers and 
Structures, Vol. 43, No. 3, 1992, pp. 565-579. 

doi: 1 0. 1 0 1 6/0045-7949(92)90289-C 

[23] Toniolo, M. D., “POST2 CFE Updates and Verification Final Report,” 
Analytical Mechanics Assoc., AMA Rept. 10-04, Hampton, VA, 
Jan. 2010. 

[24] McClinton, C. R., “X-43 Scramjet Power Breaks the Hypersonic 
Barrier, Dry den Lectureship in Research for 2006,” AIAA Paper 2006- 
0001, 2006. 

[25] Schaechter, D. B., and Levinson, D. A., “Interactive Computerized 
Symbolic Dynamics for the Dynamicist,” Journal of the Astronautical 
Sciences, Vol. 36, No. 4, Oct.-Dec. 1988, pp. 365-388. 

[26] Kane, T. R., and Levinson, D. A.., AUTOLEV 4 User’s Manual, OnLine 
Dynamics, Sunnyvale, CA, 2005. 

[27] Kane, T. R., and Levinson, D. A., Dynamics Online: Theory and 
Implementation With AUTOLEV, OnLine Dynamics, Sunnyvale, CA, 
2000. 

[28] Kane, T. R., and Levinson, D. A., Dynamics: Theory and Applications, 
McGraw-Hill, New York, 1985, pp. 158-166. 

[29] Karlgaard, C. D., Tartabini, P. V., Blanchard, R. C., Kirsch, M., and 
Toniolo, M. D., “Hyper-X Post-Flight Trajectory Reconstruction,” 
Journal of Spacecraft and Rockets, Vol. 43, No. 1, 2006, pp. 105-1 15. 
doi: 10.2514/1.12733 

[30] Tartabini, P. V., Bose, D. M., McMinn, J. D., Martin, J. G., and Strovers, 
B. K., “Hyper-X Stage Separation Trajectory Validation Studies,” 
AIAA Paper 2003-5819, Aug. 2003. 

M. Miller 
Associate Editor 


